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Abstract 

In this work we propose and analyse a numerical method for computing a family of highly 
oscillatory integrals with logarithmic singularities. For these quadrature rules we derive error 
estimates in terms of N, the number of nodes, k the rate of oscillations and a Sobolev-like regularity 
of the function. We prove that that the method is not only robust but the error even decreases, for 
■ fixed TV, as k increases. Practical issues about the implementation of the rule are also covered in this 

paper by: (a) writing down ready-to-implement algorithms; (b) analysing the numerical stability 
of the computations and (c) estimating the overall computational cost. We finish by showing some 
\ numerical experiments which illustrate the theoretical results presented in this paper. 
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1 Introduction 



in . 

This paper concerns itself with the approximation of 

CO 



2* (/) : = / /(*) M(z " «) 2 ) eMikx) dx (1.1) 



where a G [—1,1]. For the sake of simplicity, we will assume throughout this paper that k > 0, 
although the algorithm and the theoretical results can be straightforwardly adapted for k < 0. 

Our aim is to design numerical methods whose rates of convergence do not depend on k but only 
on / and the number of nodes of the quadrature rules. No information about the derivatives, which is 
very common in the approximation of oscillatory integrals (see |8] or |12| and references therein), will 
be used, which results in a simpler and less restrictive method. At first sight, a G { — 1,0, 1} could be 
the more common cases but since the analysis we develop here is actually valid for any a G [—1,1], we 
cover the general case in this paper. 

We choose in this work the Clenshaw- Curtis approach: 

ZfcX/) == f QnI(x) log((x - a) 2 ) exp(ikx) dx « l%(f), (1.2) 
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where 

Pjv 3 Q N f, s.t. (Q N f) (cos(nvr/iV)) = /(cos(nvr/iV)), n = 0, . . . , N. (1.3) 

In other words, <2at/ is the polynomial of degree iV which interpolates / at Chebyshev nodes. 

Classical, and modified, Clenshaw-Curtis rules (11 .2[) enjoy very good properties which have made 
them very popular in the scientific literature cf. [5j [201 [21] [231 ES] and been considered competitive 
with respect to Gaussian rules even for smooth integrands (we refer to [26] for an interesting discussion 
about this fact). First of all, the error of the rule is, in the worst case, like the error of the interpolating 
polynomial in L\{— 1,1). Thus, the rule is robust respect to k and it inherits the excellent approx- 
imation properties of the interpolant. On the other hand, and from a more practical view, nested 
grids can be used in the computations. Hence, if Ij. j\r(/) has been already computed, Ik,2N(f) only 
requires TV" new evaluations of /, i.e. previous calculations can be reused. Moreover, by comparing 
both approximations, a-posteriori error estimate is at our disposal almost for free. Finally, Qn/ can 
be expressed in the Chebyshev basis very fast, in about 0(N log N) operations, using FFT techniques. 

If k = 0, or if k is small enough (k < 2 has been used throughout this paper), the complex 
exponential can be incorporated to the definition of /. This leads us to consider, in the same spirit, 
the following integral and numerical approximation, 

OT) := J f(x) log((x - a) 2 ) dx « J Q N f(x) log((x - a) 2 ) dx =: X >(/). (1.4) 

This problem is also dealt with in this work since the combination of both algorithms gives rise to a 
method which can be applied to non-, mildly and highly oscillatory integrals. 

For these rule we will show that the rule converges superalgebraically for smooth functions /. 
Moreover, the error is not only not deteriorated as k increases but it even decreases as k~ l as k — > oo. 
Furthermore, for some particular values of a, which include the more common choices a £ { — 1, 0, 1}, 
the error decay faster, as A; -2 , which means that both, the absolute and relative error of the rule 
decreases )cf. Theorem I2.4|) . 

The implementation of the rule hinges on finding a way to compute, fast and accurately, the weights 

Q(k) := J T n (x)\og((x -a) 2 ) exp(ikx)dx, k> 2, (1.5) 

C ■= O0) = J l T n {x)\og{(x-af)dx (1.6) 

(T n (x) : = cos(n arccos x) is the Chebyshev polynomial of the first kind) for n = 0,1,..., N. The 
second set of coefficients (£") n is computed by using a three-term recurrence relation which we show to 
be stable. For the first set, (£"(fc)) n , the situation is more delicate. First we derive a new three-term 
linear recurrence which can be used to evaluate £"(fc). The calculations, however, turn out to be 
stable only for n < k. This could be understood, somehow, as consequence of potentially handling 
two different sources of oscillations in £"(£;). The most obvious is that coming from the complex 
exponential, which is fixed independent of n. However, when n is large, the Chebyshev polynomials, 
like the classical orthogonal polynomials, have all their roots in [—1,1]. This results in a increasing 
oscillatory behaviour of the polynomial as n — > oo. As long as the first oscillations source dominates 
the second one, i.e. as k > n, the recurrence is stable: any perturbation introduced in the computation 
is amplified very little. However, when n > k increases, such perturbations are hugely magnified, 
which makes this approach completely useless. Of course, if k is large, so should be iV to find these 
instabilities. Hence, this only causes difficulties for practical computations in the middle range, that 
is, when k is not yet very large but we need to use a relatively large number of points to evaluate the 
integral within the prescribed precision. 
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This phenomenon is not new: It has been already observed, among other examples, when computing 
the simpler integral 



(See [6] and references therein). Actually, the problem is circumvented using the same idea, the so- 
called Oliver method (cf. |18| ) which consists in rewriting appropriately the difference equation used 
before now as a tridiagonal linear system whose (unique) solution gives the sought coefficients except 
the last one which is part now of the right-hand-side. Therefore, the evaluation of this last coefficient 
has to be carried out in a different way. Thus, we make use of an asymptotic argument, namely 
the Jacobi- Anger expansion, which expresses £%-(k) as a series whose terms are a product of Bessel 
functions and integrals as in (11 . 6|> . Despite the fact that it could seem at first sight, the series can also 
be summed in about 0(N) operations. The resulting algorithm has a cost 0{N log N), cost which is 
lead by the FFT method used in the construction of the interpolant QnJ- 

Let us point out that the case of a = 0, for both the oscillatory and non-oscillatory case, has 
been previously considered in jl] using a different strategy. Roughly speaking, it relies on using the 
asymptotic Jacobi-Anger expansion for all the coefficients, no matter how large k is respect to n. 
Our approach is, in our opinion, more optimal since the algorithm is simpler to implement and the 
computational cost is smaller. 

The interest in designing efficient methods for approximating oscillatory integrals has been in- 
creased in the last years, fueled by new problems like high frequency scattering simulations cf. |4j [9j |6]. 
For instance, in the boundary integral method, the assembly of the matrix of the systems requires com- 
puting highly oscillatory integrals which are smooth except on the diagonal. Hence, after appropriate 
change of variables, we can reduce the problem to evaluate 



Typically, / is smooth except at the end-points where an integrable singularity, which could be either 
in the original integral or introduced in the change of variables, occurs. Actually, the log-singularity is 
very common since one can find it in the fundamental solutions for many differential operators in 2D, 
for instance, in the Helmholtz equation. 

Different strategies have been suggested for computing oscillatory integrals. For instance steepest 
descent methods, based on analytic continuation in the complex plane |8] or Levin methods which 
reduces the problem to solving ODE by collocation methods |15[ 119], On the other hand, we find Filon 
rules which consists in interpolating the function by a (piecewise) polynomial. Therefore, our method 
can be characterised as a Filon rule. The general case for smooth functions has been considered eg. in 
|10[ [TT| [T2| [T6| l29l [28] . Provided that the new integral with the interpolating polynomial replacing the 
original function can be computed exactly, a robust method is obtained in the sense that it converges 
as the size of the subintervals shrink to zero. Depending on the choice of the nodes we have Filon- 
Clenshaw-Curtis rules, Filon-Gaussian rules or, if the derivatives, usually at the end points, are also 
interpolated, Filon-Hermite rules. Oscillatory integrals with algebraic singularities in the integrand, 
and more general oscillators, have been considered in |17[ 113], A different approach was considered in 
[7] where the use of graded meshes toward the singularities has shown to be also efficient. Let us point 
out that this last example gives another example of the importance of having robust methods, which 
covers all possible values of k, since graded meshes can easily have very small subintervals so that the 
oscillations are reduced or even disappear. 

This paper is structured as it follows: In section 2 we derive the error estimates for the quadrature 
rule. In section 3 we deduce the algorithms to evaluate the coefficients (11. 5ft and (11. 6ft . The stability 
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of such evaluations is analysed in detail in section 4. Some numerical experiments are presented in 
section 5, demonstrating the results proven in this work. In the appendix we collect those properties 
of Chebyshev polynomials used in this paper. 



2 Error estimates for the Product Clenshaw- Curtis rule 



The aim of this section is to derive convergence estimates for the error of the quadrature rule (jl.2p - (jl.4p 
. Obviously, 

1£ N - Z£(/) = J ^ E N {x) log ((z - a) 2 ) exp(ifcr) dx 

where 

E N ■= Q N f - f, (2.1) 

and QnI is the interpolating polynomial at Chebyshev nodes cf. (jl .3j> . 

A very popular technique when working with Chebyshev polynomial approximations is to perform 
the change of variable x = cos 8. This transfers the problem to the frame of even periodic functions 
and their approximations by trigonometric polynomials. Hence, if we denote f c {9) := /(cos0) (note 
that f c is now even and 27r-periodic) we have that 

span(cosn# : n = 0, . . . N) B (Q N f) c , (Q N f) c (nir/N) = f c (mr/N), n = 0,...,N. 

Let us denote by H r (I) the classical Sobolev space of order r on an interval IcK and define 

H%:={<p€H{ oc (R)\<p = <p{2ir+.)}. 

The norm of these spaces can be characterised in terms of the Fourier coefficients of the elements as 
follows 



+ ^\n\ 2r \£(n) 



<p(n) :-- 



1 

2^ 



<p(0)exp(-m0)d0. 



If r = 0, we just have the L2(—tt,it) norm, whereas for positive integers an equivalent norm is given 
by 

1/2 

\m\ 2 de- 



(r) 



de 



The convergence estimates for the trigonometric interpolant in Sobolev spaces (see for instance 
§8.3]) can be straightforwardly adapted to prove that 



\\( E n) c \\h*{0,tt) < \\( E n) C \\H^ < C S:ro N S r ||/ c ||if£ 

where r > s > with r > tq > 1/2 and C s ^ independent of /, N and r. 
Set w(x) := (1 - x 2 ) 1 / 2 and define for (3 £ {-1, 1} 



(2.2) 



0/2 I 



i 2 (~l,l) 



|/(x)| 2 (l -x 2 f /2 dx 



1/2 



Notice in pass 



From the relations 



f(x)g(x)dx 



< ll/ll-i,t«lblli,t« 5 \\g\\i,w < V7r||5lUoo(-i,i)- 

l,w — WfcW ||/ , ||l, M = ||(/c) , || L2(0,7r) ) 



(2.3) 
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estimates (|2.2p . and the Sobolev embedding theorem |23[ Lemma 5.3.3], we can easily derive the 
following estimate: For any r > 1, 

\\EnW-i,* + N-^E'xW^ + N- l / 2 -"\\E N \\ Loo{ _ ltl) < C £ iV- r ||/ c ||^, (2.4) 

with C e depending only on e > 0. 

To prove the main result of this section we previously need some technical results we collect in the 
next three Lemmas. The first result concerns the asymptotics of 

£o (^) = J exp(ikx) log((x — a) 2 ) dx 

as k —> oo. 

Lemma 2.1 For all a € ( — 1, 1) there exists C a > so that for all k > 2 

ItfMI^C^l + togCl-o 2 ))*:- 1 . 

Moreover, for a = dbl, 

with C\ > independent of k > 2. 

Proof. The result follows from working on the explicit expression for £o(k), see fj3. 18[) - fj3. 19|) . and 
from using the limit of the functions involved as k — >• oo. We omit the proof for the sake of brevity. 

□ 

The next Lemma complements the estimates given in (12. 4ft , 
Lemma 2.2 Let En be given in (|2.ip and 

e%(x) := (2 . 5) 
x — a 

Then for all r > so > 5/2, there exists C SQ independent of f , N and r so that 

Kllw-i.i) + \\e a N \\L x (-i,i) < C S0 N s °- r \\f c \\ H r # . (2.6) 

Moreover, 



ll^llw-1,1) + \\w(e%y\\ Loo( „ lA) < C sl N s ^ r \\f c \\ H r # , (2.7) 
for r > s± > 7/2, with C S1 independent also of f , N and r. 

Proof. Recall cf. (jXjj - flSH) 

l|3^IUac(-i,i) = " 2 > ll wT nllioc(-i,i) = n > \\wTn\\ Loo (-i,i) < Cn 3 (2.8) 

where C > is independent of n. Define now, 

:= T n ( X )-T n (a) e 
x — a 

Obviously 

bs;iu 00 (-i,i)<KiUoo(-i,i)^" a - ( 2 - 9 ) 
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Besides, from (1A. 10|) (note that in the notation used there, Uj = T'- +l /{j + 1)), we derive 



n-2 



Then, using ( 12,8 



(K)'(x) = 2 Ytf + ir'T^iaK^x). 

3=0 



n-2 

\ Loo (-i,i) < 2^|(j + l)- 1 T,' +1 (a)|||^_ w || Loo( _ lil) 

3=0 
n-2 

= 2j^(i + l)(n-l-j) = -(n-l)n(n + l) < 

3=0 



y 



(2.10) 



On the other hand, for all / smooth enough 

1 

tc[n)T n , f c (n) = ■ 

71=1 



/ = / c (0) + 2^;7 c (n)r n , %{n) 



} , / c (0)exp(in0)d0 = - I f{x)T n {x)dx. (2.11) 



To prove (|2,6p . we recall first the definition of in (|2.5p . and combine (|2.1ip . (|2.9p and (|2,8p to 
obtain 



n=l 



^(-1,1) + II^IU^c-i,!) = ||X)SS( n K r , in , + II Z)SS(n)^ 

n=0 n=0 

oo 

< £|(^c(n)|(M^(-l,l) + H^IUocC-l,!)) 

1/2 



< 2 



2c 



n=l 



1/2 |- oo 

£i(S0e( 

71=1 



ra)| 2 n 5+2e 



Estimate (|2.2p proves (|2.6p . Proceeding similarly, but using the last bound in (|2.8p and (|2.10p instead, 
we prove ([277]) . □ 



Lemma 2.3 There exists C > stjc/i that for any a G [—1, 1] and for all g G H l {— 1, 1) mi/t (7(a) = 



/: 






x — a 



| 5 (s)| < Clx-al 1 / 4 ^'!!^, x€[-l,l], 
dx = Cl^'lli^. 



Proo/. Clearly, (l2~T3l) follows from (l2~T2l . 
Note first 

arcsin() — arcsin a 



C := max 
ae[-l,l] 



-a|V2 



< 00. 



Since 5(a) = 0, it follows 



</(s)ds 



< 



ds 



-,1/2 



< I arcsin x — arcsin a\ 



\g'(s)\\l-sY /2 ds 
1/2 



nl/2 



\g'(s)\ 2 w(s)ds 



<C\x-a\ l '%'\\^ 



(2.12) 
(2.13) 



The result is then proven. 



□ 



We are ready to give the main result of this section which summarises the convergence property of 
Xfc j\r(/) in terms of N, k and the regularity of /. 



Theorem 2.4 For all a G [—1, 1] there exists C a > so that for 5 G {0, 1} 

mn - KnU)\ < c a (i + k)- 5 N 8 ~ r \\f c \\ H r 



(2.14) 



for all r > 1 and k > 0. 

Furthermore, for all e > there exits C £ > such that if a = ±1 or a = anc? iV is ewen if holds 



\mf)-Zk,NU)\ ^ a(i+fc)- 2 (i+a 2 io g fc)iv 7 / 2 +— n/ c ||^ 

/or a// / G iJ^t tuift r > 7/2 + e 
Proof. Note first 



(2.15) 



-i 



En(x) log ((a; — a) 2 ) exp(i/cx) dx 



< llloe 



a 2 )\\ 1)W \\E N \\. 



■1,1V 



< CN- r \\f c \\ Hi . 



(2.16) 



where we have used (j2.4[) (see also (|2.ip ). This proves (|2.14p for S = 0. Observe that from now on we 
can assume, without loss of generality, that k > 1. 
To obtain (|2.14p for 5 = 1 we write 



mf)-ikAf)\ 



i 

< - 

- k 



+ 



(E N (x) - E N (a)) log ((x - a) 2 ) exp(ifcx) dx + E N (a)^(k) 

l 

(En(x) — Ejsf(a)) log ((x — a) 2 ) exp(i/cx) 



x=-i 



E' N (x) log ((x — a) 2 ) exp(i/cx) dx 



+2 



1 E N (x)-E N (a) 



-i 



x — a 



exp(ifcx) dx 



+ 1^(0)11^(^)1 



[R$ (k) + (fc) + R® (k) + (&)] . 



(2.17) 



For bounding the first term, we make use of (|2,12p in Lemma 12.31 with g = i^^r — Ej^(a): 
R$(k) < C||£^|| a , tt | (\h a (-l)\ + Ml)|), K(x) := |x - a) 1 / 4 log ((x - a) 2 ) < C"||^|| 1)to . (2.18) 



where we have used that h a G C [— 1, 1] for any a G [—1,1]. 
For the second term, notice that (|2.3p and (|2.4p imply 

R%\k) < || log ( - -a 2 )\U >w \\E' N \\ ltW . 
On the other hand, (|2.13p of Lemma [2.31 yields 



(2.19) 



R%\k)<C\\E' N \\ liW . 



(2.20) 



Finally, E^(±l) = and therefore R^\k) vanishes for a = dbl. Otherwise, Lemma |2 , 1 1 implies 



R%\a) < C' a \E N {a)\ < C' a \\E N \\ Loo( _ ltl) . 



(2.21) 



Bounding fl2T8l) - fl2T2TT) with (|Z3]), we derive fl2~T4l) for 5 = 1. 

To prove (j2.15p we have to perform another step of integration by parts. First, we note, that, by 
hypothesis En (a) = which implies 



R$(k) = R%\k) = 0. 



(2.22) 



Thus, we just have to estimate R N (k) and R N (k). For the first term, and proceeding as in (12 . 1 7|) . 
we obtain 



?(2)" 1 



(E' N (x) - E' N (a))\og({x - a) 2 )exp(ifca;)|^ = _ 1 



+ 

+2 

l r 

k 



-i 



E N {x) log ((x — a) 2 ) exp(ifcx) dx 



1 £^(x)-£^(a) 



l x — a 
S N \k) + S%\k) + S§\k) + S$\k) 



exp(ikx)dx +-(\E' N (a)\ \k^(k)\) 



(2.23) 



Proceeding as in (|2.19p - (|2.20p . we obtain 



S ( "\k) + S { ^(k) < C\\E' N \\ hw < VtCWwE'^i-iA) < C e N 7 ' 2 



?(3) 



Tc \\W 



where in the last inequality we have used (|2.7p of Lemma 12.21 

On the other hand, sff (k) is bounded by applying again Lemma [2. II (it is just here where the log k 
term comes up for a = dbl; note that in this case E' N (0) = for even N). 

It only remains to study S^\k). Clearly, for a = it vanishes. For a = 1, we have instead 

S N \k) < 41og2||^|| Loo( _ lil) + lim|(^(x)-^(l))log(x-l) 2 |. 

Since (l2~T2l 

\{E' N {x)-E' N {l))\og{x-lf\ < II^HUooC-i,!) |(l-^) 1/4 log(0r-l) 2 )| ->0, asx-H, 
(12. 6p implies 

S^\k) < 21og2||^|| Loo( _ lil) < CN 7 / 2 - r \\f c \\ Hi . 

The case a = — 1 is dealt with similarly. 

In brief, we have proved that for a = d=l, or for a = and ./V is even, it holds 

R^(k) < C^N 7 ^-^! + a 2 logk)\\f c \\ H r 



for e > 0, r > 7/2 + e with C e independent of N, k and /. 
Similarly, one can prove easily 



R%\k) < ~[|eft(-l)| + |eft(l)| + C||(e^|U oo( _ 1)1) 



l r, 



< C F k- 1 N 7 ^ 2+£ - r \ 



(2.24) 



(2.25) 



for all r > 7/2 + e, and C £ depending only on e > 0. 

Plugging (j2~!Mj) and (jX25)) into (j2~T7) and recalling (j2~22j) . we have completed the proof of ([2715)1 

□ 
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Remark 2.5 We will show in the last section (see Experiment 5) that the restriction of N to be even 
if a = for achieving the k~ 2 -decay of the error is really needed. In the same experiment, we can check 
that the error, specially for high values of k, is smaller for a = than for that obtained if a = 1. This 
supports empirically the fact that in the second case the log k term is certainly part of the error term 
and therefore it affects, although very slightly, the convergence of the rule. 



3 Stable computation of the weights 



When the practical implementation of the quadrature rule is considered, we face that it essentially 
reduces to find a way of evaluating £"(&) cf. fll.5|) - fjl.6|) fast and accurately. In this section we present 
the algorithms to carry out this evaluation and we leave for the next one the proofs of the results 
concerning the stability of such computations. 

For both, the oscillatory and non-oscillatory case, what we will actually compute is 



log((a; — a) 2 )U n (x) ex.p(ikx)dx 



where 



1 



1 n+l 



(3.1) 



(3.2) 



n + 1 

is the Chebyshev polynomial of the second kind and degree n. Notice that, from this definition, we 
have U-i = and, according to that, we can set 

V-i(k) :=0 

which simplifies some forthcoming expressions. From (]A.9p we have 



^(k)=^(k), C(k) = -{^(k)+^_ 2 (k)), n=l,2... 
Observe that by (1A.3[) (jA.5[) . there exists C > such that for any a E [—1, 1] and n 



\C(k)\ < 
K(k)\ < 



J \log((x-a) 2 )\dx <C, 



\U n (x)\ 2 ^/l-x 2 dx 



1/2 



log((x-a) 2 )) 



dx 



x^ 



1/2 



< c. 



(3.3) 

(3.4) 
(3.5) 



That is, these coefficients are bounded independent of n, a and k. 



3.1 The non-oscillatory case 

Recall that for k = 0, we have denoted and r/" instead of £"(0) and ??"(0) to lighten the notation. 

Assume that a/±l. Using the recurrence relation for Chebyshev polynomials cf. (jA.2|) and (|3.2p . 
we deduce for n > 1 

Vn = J U n (x) log ((x - a) 2 ) dx 

1 H 

2xU n -i(x) log((x — a) 2 ) dx — / U n -2(x) log((x — a) 2 ) dx 
-l J-i 

J (x - a)T' n {x) log((s - a) 2 ) dx + 2a<_ x - r£_ 



-2- 



(3.6) 
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Integrating by parts in the first integral, we easily see that 
-l 

(x — a)T^(x) log((x — a) 2 )dx = 



2 = 1 



X = -l 



1 fl 

2\ 



(x - a)T n {x)\og{(x - a) ) - I T n {x)\og{(x - a) ) dx - 2 / T n (x)dx 



-l J-i 



= (1 - a) log((l - a) 2 ) + (-1)" (1 + a) log((l + af) 

1 f a , if n is even, 

^[<-<_ 2 ]+ . (3.7) 

z ID, otherwise, 

where we have used that T n (l) = 1 = (— l) n T n (— 1), (13. 3ft (see also fjA.9|) ) and (lA,6p . Hence, inserting 
(|3.7p in (|3.6p and resorting appropriately the elements above, we arrive to the following two-terms 
linear recurrence 

2rVT? 71 — 1 

< = — T^-i " —T 1^-2 + it n = l,2,... (3.8) 
n + 1 n + 1 

with 

^ a . = 4 I (1 - a) log(l - a) + (1 + a) log(l + a) + ^ 2 _ - , for even n, ^ 
n + 1 \ (1 - a) log(l - a) - (1 + a) log(l + a), for odd n. 

For a = ±1, (|3.8p remains valid with 

±1 8 J log2 + , for even n, , , 

n + i I + log2, for odd n, 



which corresponds to take the limit as a — > ±1 in fj3.9|) . 
Straightforward calculations show, in addition, that 



f -( a -l)bg((a-l) 2 ) + (a + l)bg((a + l) 2 )-4, if a ^ ±1, 

\ 4 log 2 - 4, if a = ±1. ^ ' ' 



Recalling that rj" 1 = 0, we are ready to write down the first algorithm. 

Algorithm I: compute £° f° r n = 0,1, . . . , N 

1. Set r/" 1 = and compute t]q according to (|3.1ip . 

2. For n = 1,...,N 

a 2o>n n 1 a 

n + 1 n+1 

with 7° defined in (l3T9l)- (l3~T0l . 

3. Set 

£o = Vo , C = hK-V%- 2 ], n = 1 , 2, . . . , N. 

Remark 3.1 For a = the algorithm is even simpler since by parity f?2n+i = ^2n+i = and step 2 of 
the algorithm becomes 

n _ 2n-l o 



^2n — — n„ , 1 7 ?2n-2 + 



2n + l ' in ~ 2 (2n + l)(4n 2 -l)" 
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3.2 The oscillatory case 

Because of f)3.3|) . 



1 



o«(fc) + <- 2 (*0) = Q(k)= {T n {x)-T n {a))\og({x-a) 2 )ex V {ikx)dx 



-i 



+T n (a) J log((x — a) 2 ) exp(ifcr) dx. 



(3.12) 



Assume now that a/±l. Integrating by parts we derive 



x) — T n (a)) log((x — a) 2 ) exp(ikx) dx 



1 

ik 



1 

ik 



(T n (x) - T n (a)) log((x - a) 2 ) exp(ikx) 



T' n {x) log((x — a) 2 ) exp(ikx) dx — 2 



x=l 
x=-l 

1 T n (x) - T n (a) 
i x — a 



exp(ifcx) dx 



(3.13) 



(1 - T n (a)) log((l - aY) exp(^) + ((-l) n+1 + T n {a)) log((l + aY) exp(-ifc) 



ft / f n _i(x) log((x — a) 2 ) exp(i/cx) dx 



2 / 6 r n _i(x) exp(ifex) dx — 4 T n _i„j(a) / Uj(x) exp(ifex) dx 



(3.14) 



(We have applied ( lA.lOj) to write the last integral in A3. 13|> as the sum in the right-hand-side of (13. 14ft ) . 
Inserting (|3.14j) in ( 13. 12[) and using ()3.2|l we derive the following recurrence equation 



2n 



(3.15) 



where 



- [(1 - T n (a)) log((l - aY) exp(ik) + ((-l) n+1 + T n {a)) log((l + a) 2 ) exp(-ifc)] 
4 



n-2 



j=Q 



(3.16) 



with 



Pj(fc) := J Uj(x) exp(ikx)dx, j = 0, 



, n — 1. 



Let us point out that (pj(k))j =0 can be computed in 0{N) operations (see [6 
For a = ±1 we obtain the same recurrence (13. 15ft with 



4 r J log(4) exp(^i/c), if n is odd 
ifc L 1 0, otherwise 



n-2 -i „i 

2V(±1)^ ,+1 / C/ i (x)exp(iA;x)dx- / t/„_i(x) exp(ifcx)dx . (3.17) 

3=0 J - 1 J - 1 J 
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It just remains to compute ?7g (&) for setting up the algorithm. For this purpose we introduce the 
sine and cosine integral functions 

sin x f cos x — 1 

Si(i) := / dx, Ci(t) := 7 + log(t) + / dx, 

Jo x Jo x 

with 7 ~ 0.57721 the Euler-Mascheroni constant. Straightforward calculations show that 

2 



Vo t (k)=$(k) 



k L 



log(l - a 2 ) sin A; + sin(a£:) ( Ci((a + l)k)) - Ci((l - a)fc)) 
cos(afc)( Si((a + l)jfe) + Si((l - a)k)) 



2» r 
T 



,1 + a 



log (- — -) cos fe + cos(a£;) ( Ci((l - a)k) - Ci((l + a)k)) 



sm{ak) ( Si((l - a)k) + Si((l + a)k) 



(3.18) 



for a 7^ ±1, and 

vi\k) = C b tl (k) 



2i r 



(7 - Ci(2fc) + log(fc/2)) sin k - Si(2fc) cos fc 
(7 - Ci(2Jfe) + Iog(2fc)) cos k - Si(2/s) sin fc 



(3.19) 



From now on, we will denote by [x\ the floor function, i.e, the largest integer smaller than x. 

Algorithm II: computation of for n = 0, . . . , N with N < [k\ — 1 

1. Set ^(fe) = and evaluate r^(/c) according to (|3~T8]) - (|3~T9]) . 

2. Compute -yg(jfe) forn = 1, . . . , N using (l3~T6l) - (l3~T7l) . 

3. For n = 1,2,... ,7V, define 

2n 



4. Set 



r&(k)=T£(k)-—TjZ_ 1 (k)+r£_ 2 (k). 

IK 



$(k)= Vo *(k), g(k) = -[r£(k)-T&_ 2 (k)], n 



(3.20) 



1, 



,N. 



Observe that we have restricted the range for which this algorithm can be used to N < k — 1. This 
is because the recurrence relation (|3.20p . as it will be shown in the next sections, is not longer stable 
for n > k. Thus, we have to explore different ways to compute rj%(k) when n > k. 



Then assume that N > [^J — 1. Note that Algorithm II returns ?7o (&)> 



Lfcj-i- 



In order to 



compute the remaining weights we still use (|3.20p but rewriting it in a different way, namely as a 
tridiagonal system. (This is the so-called Oliver method cf. |18|). Hence, let 



A%(k) :-- 



1 



2(LfcJ+i) 

__ x 2(LfcJ+2) 1 



ik 



2N~2 
ik 



h%(k) 



'^fcj-i + Tffcj+iW 



(3.21) 
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Note that A^(k) is row dominant, which implies first that the system is uniquely solvable and next 
suggests that all the calculations become stable. This will be rigorously proven in next section. 
In short, if rj solves 

A%(k)r, = b%(k), (3.22) 

necessarily 

T 



In the definition of the right-hand-side we find r]fj(k) and ^nu^ The latter is already known. Thus, 
only the problem of finding rj^^k) remains open. For these purposes, as in [3], we will use the Jacobi- 
Anger expansion cf. [23 §2.2], (2J (9.1.44-45)]: 

oo 

exp(ifcr) = J (k) + 2 ^ i n Jn(k)T n (x), 

n=l 

where J n is the Bessel function of the first kind and order n. Hence, using (|A.7p . we derive 

ri%{k) = j Un{x) exp(ikx) log((x — a) 2 ) dx 

= Jo(fc) / U N (x)log((x-a) 2 )dx + 2 Y^i n J m (k) / U N (x)T n (x)\og({x-af)dx 

m=i J - 1 

N oo 

= Jo(fc)?7^r + ^ i m J m (k)(r]N +m — f]N-m) + i m Jm{k){ri% +m — r/m-N-2)- 

m=l m=N+l 

Observe that the coefficients rj" can be obtained from Algorithm I. By (|3.5p . in order to estimate how 
many terms are needed to evaluate this coefficient, we need to estimate how fast Ju{k) decays as 
M -> oo. We point out cf. [2j (9.1.10), (9.3.1)], 

1 (k\ M 1 / ek \M 

J-w»isi( 5 ) o^miwi) (3 - 23) 

which shows that Jjvf(fc) decreases very fast as n — )• oo. In addition, it suggests that taking ~ k terms 
in the series (|3.23p . should be enough to approximate rjN^k) within the machine precision. 
In our implementation we have taken 

TV M(k) 

V%(k)*Mk) V %(k)+J2 inj ^ k )(VN+n-ri%-n)+ £ inj n(k){v%+n-V%-N-2)- (3-24) 

n=l n =N+l 

with M{k) = 25 + [e/e/2] which has demonstrated to be sufficient for our purposes. 

Algorithm III: compute for n = [k\ , . . . , N 

1. Construct h^(k) using 

(a) ^nu_ returned in Algorithm II 

(b) 7£(Jfe) for n = [A;J + 1, . . . , A defined in (|3TT5]l - ([3"37]l . 

(c) r)x(k) evaluated with the sum f|3. 24|) . 
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2. Construct the tridiagonal matrix A^(k) denned in (|3.2ip and solve 

A a N (k)v = h a N (k). 

Set 

vw-i+M = (v)t, e = i,...,N-ik\. 

3. Set 

C(k) = ~[<(A0-<- 2 (A0L n=[k\,...,N. 
On the computational cost 

Certainly, one could use (|3.24p for computing all the coefficients (ry°(/c)) n , as it was proposed in [3] 
(for a = 0). However, this choice results in a more expensive algorithm. By restricting this approach 
to the last coefficient, and only if N > k, we can speed up the algorithm since all the terms but the 
last one, are computed by solving a tridiagonal system which can be done in 0{N — \ k\) operations 
by Thomas algorithm. 

The vector (7"(fc)) n= i can be also constructed very fast. Hence, note that the bulk part in (|3.9jl is 
the convolution of the vectors 

which can be done in 0(Nlog(N)) operations by using FFT. (For a G { — 1,0, 1} this could be achieved 
even faster from (|3.9p . since T n (±l) = (±l) n and T n (0) = 1 if n is even and otherwise). 

Another possible bottleneck of the algorithm could be found in the evaluation of the Bessel functions 
J n (k). Let us show how it can be overcome. We recall that the Bessel functions obey the recurrence 
relation 

2n 

Jn+i(k) - —J n {k) + Jn-l(k) = 0. (3.25) 

Notice in pass that it is very similar to that obtained in (|3.15j) for evaluating our coefficients. Thus, 
we can exploit these similarities to get a faster evaluation of these functions: Once Jo(k) and J\{k) 
are evaluated by usual methods, (13. 25ft can be safely used for evaluating J n (k) for n < [k\. For the 
remainder values, i.e. for n > [k\ + 1, we make use of the Oliver approach and solve 

' 2(Lfcj+i) x 

J* 2(LfcJ+2) 1 

k 

1 2M(k) 
1 k . 

The asymptotics (|3.23p can be used to approximate JM(k)+i(k), which gives even better results that 
setting simply JM(k)+i(k) ~ 0. The evaluation turns out to be stable just for the same reasons that 
ensure the stability of Algorithms II and III (see next section). 

4 Numerical stability 

We analyse the stability of the algorithms separately in three propositions and collect the stability 
results for Algorithms II and III, when they work together, in a theorem which ends this section. 

The (usually small) parameter Sj > will be used in this section to represent any possible perturba- 
tion occurring in the evaluation such as round-off errors or errors coming from previous computations. 



J lk\~ 


n(*) 




-J[k\-i(k) 


J[k}- 


1-2 W 







Jm( 






-■h4(k)+i(k) 
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Theorem 4.1 Let := (eoj £ 1i ■ ■ ■ i £ N) £ with \\s\\oo < e an d define the sequence 

V-i = V-i = °> Vo =Vo + £ o, 

C = 7n + " + £ n, for n = 1,2, ... . 



TTien /or all N > 



I Q a,£| ^- 

WjV — I < 



1 N 

j=0 



e < \{N + 2){N + 3)e. 



Proof. Clearly, 



where On is given by 



N 



3=0 



:= 0, <5f: 



(i) ._ 



-3 ) u n 



2an 



-S 



U) 



n-\ 



n + l n ~ L n + l 

It is easy to check, using (1A.2|) . that the solution of the problem above is given by 



n = j + l,j+2,... (4.1) 



Therefore, using (IA.3|) 



l a a,£\ 



< 



1 



A? 



3=0 



N + l 
U N + 2)(N + 3)e. 



£ < 



1 



N + l 



N 



Y^(j + l)(N-j + l) 

3=0 



The proof is now finished. 



□ 



Remark 4.2 In view of this result, we conclude that theoretically a = turns out to be the most stable 
case. Indeed, since U2j(0) = (— l) 3 , 

\V2N - r] 2N \ < + I) - m ■ 

3=0 

(Note that V2j+i = ® an< ^ therefore only (r/^Oj has to be considered). 

On the other hand, a = +1 are precisely the most unstable cases, since \Uj(+l)\ = j + 1. We point 
out, however, that in practical computation the algorithm has demonstrated, see section 5, that: (a) the 
computation is stable for a E [—1,1], much better than that theory predicts; (b) the error observed for 
a = is a little smaller than that for a = +1. 

Next we consider the stability of Algorithms II and III, i.e., the oscillatory case. 
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Proposition 4.3 Let N < k — 1 and set e = (e , . . . , £jv) G C^ 44 wii/i H^Hoo < £ ari ^ consider the 
sequence 



v-m = 

Then, for all < n < k — 1. 



7nW " ^-l(*0 + <-2 W + n = 1, 2, 



! 4 (n+l)^ 1 / 2 



Therefore, for all n < k — 2. 



V r(k)-Vn(k)\< 



1 + 



3 (fc 2 - (n + l) 2 ) 1 /* 

2 3/4 



+ If/'k 1 



/2 



e. 



(4.2) 

(4.3) 
(4.4) 



whereas for k — 2 < n < k — 1, i.e., for n = [k\ — 1, 

l^-iW-^j-i( fc )l^[ 4 + 27/4fc5/4 ] £ - 
Proof. As before, it suffices to study the sequence 

5-i = 0, 6q = eo 
2n 

o~ n = — rp^n-i + 5„-2 + e n , n = l,2, ... 

IK 

We refer now to j6j Th. 5.1] where the stability of this sequence is analysed and whose proof can be 
straightforwardly adapted to derive (|4.2p 

To prove (14. 3ft . we observe that ( 14. 2ft implies that for n < k — 2, 



AT. 



|<' £ (fc)-7£(fc)| < 



1 + 



n 



+ 1)^/2 



3((n + 2) 2 -(n + l) 2 )V4 

97/4 



e < 



1 + 



4 (n + l)*: 1 / 2 
3 (2n + 3)V4 



If n = \k\ — 1, we can use (|4.3p as follows 

\vr(k)~vm\ < £ + T K^(k) - V%-i(k)\ + Kf 2 (k) - vum 

< [4 + 2 • ^n 3 ^ 2 + ^(n - l) 3 / 4 fcV2] e < [ 4 + 2 7 / 4 n 3 / 4 A: 1 / 2 ] e 

< [4 + 2 ? / 4 A: 5 / 4 ] e . 
Bound (|4.4p is now proven. 



□ 



The stability of Algorithm III is consequence of the next result. 
Proposition 4.4 Let N > k and consider the solutions of the original and perturbed systems 

A%(k)r 1 = h%{k), (A%{k) + AAfjikM = bft(fe) + Ah a N (k). 
Then, if (k + 2)\\AA%(k)\\ oa <2, it holds 

k + 2 



\rf - m\oo < 



+ 2)||AA° 
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Proof. A classical result in stability theory for systems of linear equations (see for instance O Th. 
8.4]) states that 



\v -v\ 



< 



l-\\AA-(k)\\ 00 \\{A%(k))- 1 \\ oo 
Thus, we just have to estimate || (^y-(fc)) 1 || o- Let 

\k\ + 1 

D N (k) 



\Ab%(k)\\ 00 + ||A^(fc)||ooN|U • (4.5) 



2 

ik 



\k\ +2 



N — 1 



K N (k) 



ik 



ik 



ik 



2(LfeJ+i) 



2(^+2) 

mm 

ik n ik 

~2(LfcJ+2) U 2(LfcJ+4) 



•ik 







2(JV-2) 

Let 7 denote the identity matrix. Clearly, it holds 

A a N {k) = {I + K N {k))D N {k) => (A^(fc))" 1 = (7J JV (A:))" 1 (7 + 7C i v(A:))- 1 . 

Notice also 

fc + 1 



+ 



1 fc 

< r + 



2([fcJ+l) 2{[k\+2,) 2 2(/c + 2) k + 2' 
Collecting these inequalities, we conclude 

I,, ^ WiDNik))- 1 ^ k + 2 



\\[AMk)y L \\ < i 

Inserting (|4.6D in (|4.5jl the result is proven. 



ll^(*0llc 



(4.6) 
□ 



The perturbation AA^(k) in the matrix is essentially round-off errors. Since A^(k) is a tridiagonal 
matrix we can safely expect (k + 2) [| (A;) || oo << 1. 

The last result of this section states the numerical stability of the algorithm in the oscillatory case 
and is result of combining appropriately propositions 14.31 and 14.41 



Theorem 4.5 With the notations of Propositions \J73\\4-4 it holds 
(a) For any v 6 (0, 1) there exists C v , depending only on v, so that for N < vk, 

max AT K' e (k)-rj%(k)\ <C v Ns. 

n=0,...,N 



(4.7) 



with C v independent of k and N . 
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«&» 




n = 


10 


1.11E-16 


4.33E- 


16 


n = 


200 


1.11E-16 


1.07E- 


15 


n = 


400 


1.11E-16 


1.31E- 


15 



Table 1: Absolute and relative errors in £° for different values of n 





e abs( 


n) 




n = 


10 


5.83E 


-16 


3.70E- 


14 


n = 


200 


5.83E 


-16 


7.92E- 


14 


n = 


400 


5.83E 


-16 


3.90E- 


13 



Table 2: Absolute and relative errors for different values of n 

(b) There exists C > independent of k and a so that 

max K> e (A:)-f£(A0| < Cfc 5 / 4 e. (4.8) 
n=Q,...,|_fcJ— 1 

(c) For N > k, if the following conditions are, in addition, satisfied 

UA^AOIU < e, HAb^^lU < e + \^_ x [k) - V^Wl 
and e < (k + 1) — , i/ten it holds 

max |r£' e (fc)-77?(fc)| < (k + 2)(l + l^j.^fe) - V^-M + \\v\\oo)e (4.9) 

< C'k g / 4 e (4.10) 

where C is independent of k and N . 

Proof. Estimates (j4.7p ~ (j4.8p follow from Proposition 14.31 For item (c) (|4,9p is a consequence of 
Proposition 14.41 Finally, (|4.10p is obtained by applying (|4.8p . which bounds the last term in (|4.9p . and 
using that ||//||oo is uniformly bounded independent of and a cf. (|3.5p . □ 



Let us emphasise that, in our computations, (I4,10p has been demonstrated to be very pessimistic. 

5 Numerical Experiments 

We collect in this section some numerical experiments to illustrate the theoretical results presented in 
this paper. The implementation of the rule, for a = 0, —1 is available in [T]. 



5.1 Stability for £ 



We have computed here £™ for n = 0, . . . , 100 using the implementation of our method in Matlab. Next, 
we compare the numerical results with that obtained using symbolic calculations in Mathematica, which 
will be denoted by ymb . The evaluation of these expressions is done using (very) high arithmetic 
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N\k 


10 


20 


40 


80 


160 


1 


1.39E- 


17 


1.04E- 


17 


1.30E- 


18 


4.34E-19 


2.71E- 


320 


10 


1.33E- 


15 


2.22E- 


16 


2.78E- 


17 


2.78E-17 


0.00 




on 
ZU 


1.67E- 


16 


6.66E- 


16 


0.00 




z. IB X 1U 


6.94E- 


18 


40 


2.78E- 


17 


1.39E- 


16 


1.11E- 


15 


4.16E-17 


0.00 




80 


2.78E- 


17 


1.39E- 


17 


5.55E- 


17 


1.11E-15 


2.08E- 


17 


160 


0.00 




1.39E- 


17 


2.78E- 


17 


7.63 x 10~ 17 


1.55E- 


15 


N\k 


10 




20 




40 




80 


160 




1 


1.58E- 


-16 


1.64E- 


-15 


6.18E- 


-16 


2.70E-16 


1.28E- 


16 


10 


6.48E- 


-16 


4.63E- 


-16 


1.63E- 


-16 


3.42E-16 


0.00 




20 


3.91E- 


-16 


3.61E- 


-16 


0.00 


3.14 x 10 


1.76E- 


16 


40 


1.65E- 


-16 


6.60E- 


-16 


6.76E- 


-16 


3.56E-16 


0.00 




80 


3.48E- 


-16 


1.66E- 


-16 


5.30E- 


-16 


7.66E-16 


8.18E- 


16 


160 


0.00 


3.48E- 


-16 


6.63E- 


-16 


1.46 x 10~ 15 


1.23E- 


15 



Table 3: Absolute (top) and relative (below) error in computing £°(fc) 



precision to keep the round-off errors well below the significant digits returned in our implementation 
in Matlab. 
We present 

|£« _ £ a > s y mb | 

e« bs( A0 = n m^ |e° - £f ymb |, eUN) = ^ J ^ mb| 

in Tables [1] (for a = 0) and [2] (for a = 1) for different values of n. We clearly see that for all n the 
error is very close to the machine's unit round off and that the results returned for a = are slightly 
better than that for a = 1. This should indicate that Theorem 14. II is sharp (see also Remark I4.2|) . 

5.2 Stability for 

As before, we compare here the values of rj%(k) computed by our code with that returned by Mathe- 
matica. The results are shown for a = in Table [3] and for a = 1 in Table [4] 

It is worth mentioning that in our implementation in Matlab we face an annoying bug. Algorithm 
II (and therefore indirectly Algorithm III) makes use of the sine and cosine integral functions (Si and 
Ci in our notation) just for evaluating the first coefficient rj^k). These functions are only included in 
Matlab as part of the symbolic toolbox, and therefore it is not presented in all distributions. Moreover, 
any call to these functions consumes a significant CPU time because of the own nature of the symbolic 
toolbox. Hence, in some of our experiments we observed that when using the built-in functions almost 
half of the CPU time was consumed in performing these two evaluations. 

Thus we wrote our own implementation for sine and cosine integral functions. The evaluation is 
accomplished by a combination of asymptotic expansion for large arguments [21 (5.2.34)-(5.2.35)] and 
a sum of Bessel functions of fractional order for small and moderate arguments [2J (5.2.15)]. Despite 
our efforts, our code introduces a very small error in the last or in the last but one significant digits. 
However, such errors only affect the first few coefficients very slightly and do not propagates to the 
rest of coefficients. Hence, it gives us a unwanted proof of the stability of the algorithm. 
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N\k 


10 


20 


40 


80 


160 


1 


3.86E- 


16 


1.39E- 


17 


2.95E- 


16 


2.78E- 


17 


1.39E- 


17 


10 


2.54E- 


15 


2.24E- 


16 


4.79E- 


16 


2.86E- 


17 


1.55E- 


17 


20 


2.22E- 


17 


1.56E- 


15 


1.25E- 


15 


2.08E- 


17 


4.39E- 


17 


40 


1.03E- 


16 


3.71E- 


17 


4.10E- 


15 


1.67E- 


16 


1.12E- 


16 


80 


1.81E- 


17 


6.35E- 


17 


2.70E- 


17 


1.60E- 


15 


1.31E- 


16 


160 


1.32E- 


17 


1.86E- 


17 


1.00E- 


16 


5.02E- 


16 


8.68E- 


16 


N\k 


10 




20 




40 




80 




160 




1 


5.47E- 


16 


3.05E- 


17 


1.35E- 


15 


1.89E- 


16 


2.15E- 


16 


10 


4.55E- 


15 


5.19E- 


16 


2.00E- 


15 


1.87E- 


16 


1.95E- 


16 


20 


7.89E- 


16 


4.42E- 


15 


3.57E- 


15 


1.04E- 


16 


4.97E- 


16 


40 


1.18E- 


14 


3.84E- 


15 


1.89E- 


14 


4.12E- 


16 


5.95E- 


16 


80 


6.94E- 


15 


2.23E- 


14 


9.93E- 


15 


1.21E- 


14 


6.56E- 


16 


160 


1.73E- 


14 


2.28E- 


14 


1.27E- 


13 


6.02E- 


13 


1.01E- 


14 



Table 4: Absolute (top) and relative (below) error in computing 

5.3 Experiments for an oscillatory integral 

Let i 

I a (k) := f log [(x — a) 2 ) exp(ikx) dx (5.1) 

J_i x + x + 1 

We have computed the errors returned by our numerical method for different values of k, N and for 
a = (Table [5]) and a = 1 ( Table |6|). As exact value we just have used that returned by the rule when 
a huge number of points is used. 

Several facts can be observed right from the beginning. First, the convergence is very fast: with 
modest values of we get approximations with an error of the same order as the round-off unity. 
Second, if we read the Table by rows, we clearly see that in both cases the error decrand it probably 
deserves more research on this topic eases as k~ 2 only for even N, whereas for odd N and a = the 
error decreases only as (and therefore, the relative error keeps bounded independent of k). 

Such phenomenon doest not occur when a = 1, i.e, when the logarithmic singularity occurs at the 
end of the interval (See Table [6]). 

5.4 Non-smooth functions 

In this last experiment we run our code to compute 

h(k,a,l3) := J (1 + x) 13 log ({x - a) 2 ) exp(ifcx) dx, (5.2a) 

I (k,a,P) := J |l/2 + x| /3 log((x-a) 2 )exp(i/cx)dx, (5.2b) 

for a G { — 1,0} and /3 G {1/2,3/2} to analyse how precise are the regularity assumptions in the 
hypothesis Theorem 12.41 

We expect the convergence of the rule to be faster for the first integral since, after performing the 
cosine change of variables, |1 + cos Of G Hf +1/2 6 whereas 1 1/2 + cos6\' 3 G H^ +1/2 e . The regularity 
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N\k 





10 


10 2 


10 3 


10 4 


10 5 


11 


1.71E- 


-03 


4.00E- 


-03 


1.75E- 


04 


1.82E- 


05 


1.83E- 


06 


1.83E- 


-07 


12 


4.56E- 


-05 


3.28E- 


-04 


1.44E- 


06 


1.37E- 


08 


1.37E- 


10 


1.37E- 


-12 


23 


1.65E- 


-08 


2.56E- 


-08 


4.80E- 


09 


3.89E- 


10 


3.80E- 


11 


3.80E- 


-12 


24 


2.96E- 


-10 


8.24E- 


-09 


9.93E- 


10 


9.09E- 


12 


9.09E- 


14 


9.08E- 


-16 


47 


6.66E- 


-16 


1.11E- 


16 


8.97E- 


17 


1.29E- 


17 


1.08E- 


19 


1.36E- 


020 


48 


6.66E- 


-16 


2.73E- 


16 


8.85E- 


17 


1.26E- 


17 


1.08E- 


19 


2.71E- 


020 


N\k 







10 




10 2 




10 3 




10 4 




10 5 




11 


9.38E- 


-04 


5.49E- 


-03 


2.78E- 


-03 


2.89E- 


-03 


2.91E- 


-03 


2.91E- 


-03 


12 


2.50E- 


-05 


4.50E- 


-04 


2.28E- 


-05 


2.18E- 


-06 


2.17E- 


-07 


2.18E- 


-08 


23 


9.04E- 


-09 


3.51E- 


-08 


7.61E- 


-08 


6.20E- 


-08 


6.05E- 


-08 


6.04E- 


-08 


24 


1.63E- 


-10 


1.13E- 


-08 


1.58E- 


-08 


1.45E- 


-09 


1.45E- 


-10 


1.45E- 


-11 


47 


3.66E- 


-16 


1.52E- 


-16 


1.42E- 


-15 


2.05E- 


-15 


1.73E- 


-16 


2.17E- 


-16 


48 


3.66E- 


-16 


3.75E- 


-16 


1.40E- 


-15 


2.01E- 


-15 


1.73E- 


-16 


4.32E- 


-16 



Table 5: Absolute (top) and relative (bottom) errors for integral (|5.ip with a = 



N\k 





10 


10 2 


10 3 


10 4 


10 5 


11 


1.81E- 


-05 


8.89E- 


-04 


3.04E- 


-05 


5.04E- 


-07 


6.33E- 


-09 


7.90E- 


11 


12 


2.43E- 


-06 


7.72E- 


-05 


8.94E- 


-06 


1.74E- 


-07 


1.77E- 


-09 


2.15E- 


-11 


13 


4.21E- 


-11 


2.60E- 


-11 


1.50E- 


-09 


5.51E- 


-12 


1.25E- 


-13 


1.48E- 


-15 


24 


5.25E- 


-11 


4.91E- 


-11 


1.89E- 


-09 


1.84E- 


11 


2.81E- 


-13 


3.55E- 


15 


47 


1.04E- 


-18 


7.85E- 


-17 


9.22E- 


-17 


2.47E- 


-17 


2.09E- 


-18 


1.10E- 


-19 


48 


7.31E- 


-17 


8.89E- 


-17 


9.17E- 


-17 


2.17E- 


-17 


1.89E- 


-18 


1.12E- 


-19 


N\k 







10 




10 2 




10 3 




10 4 




10 5 




11 


8.09E- 


-04 


3.07E- 


-03 


1.00E- 


-03 


1.71E- 


-04 


1.26E- 


-05 


1.27E- 


■06 


12 


1.09E- 


-04 


2.67E- 


-04 


2.94E- 


-04 


5.92E- 


-05 


3.53E- 


-06 


3.45E- 


-07 


23 


1.89E- 


-09 


9.00E- 


-11 


4.92E- 


-08 


1.87E- 


-09 


2.49E- 


-10 


2.39E- 


11 


24 


2.30E- 


-09 


1.70E- 


-10 


6.22E- 


-08 


6.26E- 


-09 


5.61E- 


-10 


5.72E- 


11 


47 


4.67E- 


-17 


2.71E- 


-16 


3.03E- 


-15 


8.38E- 


-15 


4.17E- 


-15 


1.77E- 


15 


48 


3.27E- 


-15 


3.07E- 


-16 


3.02E- 


-15 


7.37E- 


-15 


3.78E- 


-15 


1.80E- 


-15 



Table 6: Absolute (top) and relative (bottom) errors for integral (15. ip with a = 1 
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of the transformed function is precisely what appears in the estimate of Theorem 12,41 (function f c in 
the right-hand-side). 

We show in Tables [TlfTOl the error of the rule for different values of k and N. (The exact integral 
was computed by using the Clenshaw-Curtis rule on graded meshes towards the singular points, cf. 
[14j). Clearly the errors are in almost all cases smaller for (|5.2ap than for (I5.2bp , 

It is difficult to estimate the order of convergence of the rule because it becomes chaotic as k 
increases in such a way that the larger is k, the bigger has to be ./V to make the error decay steady to 
zero. This is specially noticeable in Tables 151 and 1101 One can observe, however, that the error in the 
first columns of [7] decreases as N~ 3 ~ 2 ^, in Table [9] the rule converges approximately as _/V~ 2 - 5 ~^ and 
as N- 1 -? for Tables MB 

If we read the table by rows, we can detect that the 0(k~ 2 ) decay of the error occurs only in Table 
[7]and in[9]for (3 = 3/2. Only for the first integral fj5. 2a|) with /3 = 3/2, this property has been rigorously 
proved since in the notation of Theorem l2.4l f r . £ H^J 2 £ . There is however no theoretically justification 
for the other cases and it certainly deserves more attention to study if the regularity assumptions can 
be relaxed for a = 0. 

On the other hand, the error does not behave as 0(k~ 2 ) in Table [8] although for /3 = 3/2 Theorem 
12.41 should imply such decay of the error. We think that the very irregular convergence of the rule in 
this case could force N and k to be larger to observe it. 



A Some relevant properties for Chebyshev polynomials 

For the sake of completeness we present in this section those properties of Chebyshev polynomials we 
have used in this work. These results can be found in many classical text books on special functions 
or Chebyshev polynomials (see for instance |2| Ch. 22] or |22|). 

From the definitions of the Chebyshev polynomials of first and second kind we have the relations 

TJcose) = cosn6, U n (cos9) = — 1— r'(cosfl) = sin ( n + ^ (A.l) 

n + 1 sin 6 

As a byproduct, one can deduce that if n is even (respectively odd), so are T n and U n . Note that 
as usual in this work, we have taken U-\ = 0, which is also consistent with fjA. 1 1) . Both families of 
polynomials obey the recurrence relation 

P n+ i(x) = 2xP n (x) - P n _i(x) (A.2) 

but with, obviously, different starting values, simply 7q(x) = 1, 7\{x) = x and U-\{x) = 0, Uq(x) = 1 
respectively. 

From (jA.ip we easily deduce 

m ; / „n sinn# d / sin n9\ 

Z£ cos0 = n— — , sin 07% cos 6) = n —[—-\. 

sin at) \ sin 6 J 

Therefore, 

11^11x^(0,1) < T n(l) = 1, 11^11^(0,1) = ^pYll T n+lllLoo(-l,l) = 71 + 1 ( A ' 3 ) 

and (recall that w(x) = y/l — x 2 ) 

ll*< +1 || Loo(0 ,i) = (n + 1), 11^111^(0,1) < C(n + If (A.4) 
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N\k 





1 


10 


10 2 


10 3 


10 4 


10 5 


11 


2.59E- 


03 


2.59E- 


03 


6.90E- 


03 


2.80E- 


04 


2.81E- 


05 


2.83E- 


06 


2.83E- 


07 


12 


1. I ZEj — 


C\A 

U4 


1. 1 Zhi — 


C\A 

U4 


A QOT7 1 
4.0ZH/ — 


AO 

Uo 


o.oorL/ — 


Uo 


Q OQT7 


U / 


z.y (tj — 


uy 


O HIT? 


"1 "1 
11 




2.86E - 


04 


2.86E- 


04 


2.96E- 


04 


4.80E - 


05 


6.48E - 


06 


6.56E- 


07 


6.57E- 


08 


24 


1.07E- 


05 


1.07E- 


05 


1.29E- 


04 


2.24E - 


05 


1.90E- 


07 


1.53E- 


09 


1.50E- 


11 


47 


3.36E - 


05 


3.36E- 


05 


3.18E- 


05 


3.04E - 


05 


1.50E- 


06 


1.57E- 


07 


1.58E- 


08 


48 


6.64E - 


07 


6.64E - 


07 


6.91E- 


06 


1.20E- 


05 


1.05E- 


07 


7.61E- 


10 


7.78E - 


12 


95 


4.07E - 


06 


4.07E - 


06 


3.91E- 


06 


5.62E- 


05 


4.14E- 


07 


3.84E - 


08 


3.86E- 


09 


96 


4.13E- 


08 


4.13E- 


08 


4.17E- 


07 


9.36E- 


05 


5.52E- 


08 


3.03E - 


10 


4.05E - 


12 


N\k 







1 




10 




10 2 




10 3 




10 4 




10 5 




11 


5.07E - 


05 


5.05E- 


05 


1.34E- 


04 


5.33E - 


06 


5.51E- 


07 


5.56E- 


08 


5.56E- 


09 


12 


2.66E — 


06 


2.65E — 


06 


8.00E — 


05 


4.74E — 


07 


4.72E — 


09 


4.78E — 


11 


4.76E — 


13 


93 
zo 


1.32E- 


06 


1.32E- 


06 


1.37E- 


06 


2.34E- 


07 


2.99E- 


08 


3.03E - 


09 


3.03E - 


10 


24 


4.41E- 


08 


4.41E- 


08 


5.47E - 


07 


9.01E- 


08 


6.59E- 


10 


6.27E- 


12 


6.17E- 


14 


47 


3.74E - 


08 


3.74E - 


08 


3.55E- 


08 


3.35E- 


08 


1.69E- 


09 


1.76E- 


10 


1.76E- 


11 


48 


7.01E- 


10 


7.01E- 


10 


7.38E - 


09 


1.25E- 


08 


1.01E- 


10 


8.11E- 


13 


7.83E - 


15 


95 


1.11E - 


09 


1.11E- 


09 


1.07E- 


09 


1.54E - 


08 


1.12E- 


10 


1.05E- 


11 


1.06E- 


12 


96 


1.10E - 


11 


1.10E- 


11 


1.12E- 


10 


2.50E- 


08 


1.36E- 


11 


9.53E- 


14 


9.99E- 


16 



Table 7: Errors of the quadrature rule for integral fj5.2a() with a = 0, /3 = 1/2 (top) and f3 = 
(bottom) 



N\k 





1 


10 


10 2 


10 3 


10 4 


10 5 


11 


4.02E - 


03 


4.02E - 


03 


1.86E- 


03 


4.99E - 


03 


3.34E- 


04 


1.56E- 


05 


6.38E- 


07 


12 


3.15E- 


03 


3.15E- 


03 


3.54E - 


03 


4.73E - 


03 


3.29E- 


04 


1.56E- 


05 


6.37E- 


07 


23 


5.11E- 


04 


5.11E- 


04 


5.11E- 


04 


2.17E- 


03 


2.79E- 


04 


1.48E- 


05 


6.27E- 


07 


24 


4.54E - 


04 


4.54E - 


04 


4.53E - 


04 


1.96E- 


03 


2.74E - 


04 


1.47E- 


05 


6.26E- 


07 


47 


6.83E - 


05 


6.83E- 


05 


6.83E - 


05 


2.20E - 


04 


1.88E- 


04 


1.32E- 


05 


6.04E - 


07 


48 


6.43E - 


05 


6.43E - 


05 


6.43E - 


05 


1.90E- 


04 


1.85E- 


04 


1.31E- 


05 


6.03E - 


07 


95 


9.28E- 


06 


9.28E- 


06 


9.28E- 


06 


1.02E- 


05 


4.78E- 


05 


1.04E - 


05 


5.61E - 


07 


96 


9.00E - 


06 


9.00E- 


06 


9.00E- 


06 


1.38E- 


05 


4.59E- 


05 


1.04E - 


05 


5.61E - 


07 


7V\fc 







1 




10 




10 2 




10 3 




10 4 




10 5 




11 


8.62E - 


06 


8.66E- 


06 


4.05E- 


05 


3.66E- 


05 


1.02E- 


06 


1.65E- 


08 


2.20E- 


10 


12 


5.19E- 


06 


5.21E- 


06 


2.07E - 


05 


3.09E- 


05 


9.05E- 


07 


1.50E- 


08 


2.01E- 


10 


23 


1.31E- 


07 


1.31E- 


07 


1.34E- 


07 


6.06E- 


06 


3.56E- 


07 


7.13E- 


09 


1.02E- 


10 


24 


1.01E - 


07 


1.01E- 


07 


1.04E- 


07 


5.30E - 


06 


3.32E- 


07 


6.78E- 


09 


9.74E - 


11 


47 


1.52E- 


09 


1.52E- 


09 


1.55E- 


09 


2.96E- 


07 


9.38E- 


08 


2.87E- 


09 


4.68E - 


11 


48 


1.30E- 


09 


1.30E- 


09 


1.33E- 


09 


2.65E- 


07 


8.95E- 


08 


2.78E- 


09 


4.57E- 


11 


95 


1.64E - 


11 


1.64E - 


11 


1.67E- 


11 


2.68E- 


09 


1.23E- 


08 


9.59E- 


10 


2.04E - 


11 


96 


1.65E- 


11 


1.65E- 


11 


1.67E- 


11 


1.76E- 


09 


1.18E- 


08 


9.41E- 


10 


2.02E- 


11 



Table 8: Errors of the quadrature rule for integral (15.2ap with a = —1, (3 = 1/2 (top) and /3 = 
(bottom) 
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N\k 





1 


10 


10 2 


10 3 


10 4 


10 5 


11 


2.13E- 


02 


2.15E- 


02 


540E- 


02 


2.96E- 


03 


8.24E - 


05 


1.31E- 


05 


1.27E- 


06 


12 


^ 71 TT 


no 
uz 


0. I Zjej — 


no 
uz 


1 01 TT 1 
1.Z1H/ — 




I.DO-Cj — 


Uo 




uo 


1.0 a Hi — 


Oft 

uu 


D.4DrL — 


Uo 


23 


5.85E - 


03 


5.86E - 


03 


6.49E - 


03 


240E- 


03 


2.81E- 


05 


4.38E - 


06 


3.79E- 


07 


24 


Z.15ii — 


02 


o i rp 


no 
02 


O A flD 


no 
Uz 


l.DDhj — 


no 

Uo 


5.o6rj — 


n ct 
05 


1 71 T7 1 

1. f lrL/ — 


06 


5.46hj — 


no 

Uo 


47 


l.ooJl/ — 


no 

(Jo 


l.odih — 


no 

Uo 


i. (ohj — 


no 
Uo 


1 

l.Of -CJ — 


no 
Uo 


4.0lHi — 


Uo 


Z.Z4H/ — 


Uo 




U / 


48 


7 '70TT' 

I . ( oih — 


no 

Uo 


/ . ( OCj — 


no 

Uo 


O.Ul-Cv — 


Uo 


1 Cfll? 

l.oyrlj — 


no 
Uo 


i5.o9±L — 


nK 
Uo 


1. ( zHj — 


n£ 
Uo 


0.4 / Ji — 


no 

Uo 


95 


d.11±!j — 


C\A 

U4 


d.11±!j — 


C\A 

U4 


r OOT7 
D.VZEj — 


C\A 

U4 


1 flCT7 


no 
Uo 


5.1 /rj — 


Uo 


1 Q .1 T? 
1.o4Hj — 


n£ 
Uo 


C? (?rri 

D.DOrj — 


no 

Uo 


96 


Z. /OHj — 


no 

Uo 


Z. 1 oHj — 


no 
Uo 


O HIT? 
Z. 1 (Ej — 


no 
Uo 


b.oohj — 


no 
Uo 


r A OT7 

0.4Z1L — 


nK 
Uo 


1. f ZH/ — 


n£ 
Uo 


r A Q T7 1 

0.4oH/ — 


Uo 


N\k 







1 




10 




10 2 




10 3 




10 4 




10 5 




11 


2.08E - 


03 


2.14E- 


03 


6.62E- 


03 


2.23E- 


04 


2.Q1E- 


05 


2.02E- 


06 


2.02E- 


07 


12 


1.39E- 


03 


1.40E- 


03 


5.83E- 


03 


3.02E- 


05 


2.47E - 


07 


2.47E - 


09 


2.65E- 


11 


23 


2.20E- 


04 


2.21E- 


04 


2.96E- 


04 


5.41E- 


05 


2.79E- 


06 


2.88E- 


07 


2.88E- 


08 


24 


2.87E - 


04 


2.87E- 


04 


4.17E- 


04 


2.47E - 


05 


1.02E- 


07 


8.26E- 


10 


1.00E- 


11 


47 


2.88E- 


05 


2.89E- 


05 


2.91E- 


05 


3.11E- 


05 


3.53E- 


07 


4.29E- 


08 


4.28E- 


09 


48 


5.32E - 


05 


5.32E- 


05 


5.87E- 


05 


2.34E - 


05 


7.33E- 


08 


1.45E - 


10 


3.15E- 


12 


95 


4.40E - 


06 


4.40E - 


06 


4.17E- 


06 


2.90E- 


05 


3.98E- 


08 


6.83E- 


09 


6.79E- 


10 


96 


9.51E- 


06 


9.51E- 


06 


9.76E - 


06 


5.08E- 


05 


7.70E- 


08 


1.22E- 


10 


5.94E - 


13 



Table 9: Errors of the quadrature rule for integral (I5.2bp with a = 0, f3 = 1/2 (top) and /3 = 
(bottom) 



N\k 





1 


10 


10 2 


10 3 


10 4 


10 5 


11 


1.74E - 


02 


1.75E- 


02 


2.53E- 


02 


1.63E- 


03 


5.63E- 


05 


1.73E- 


06 


5.49E - 


08 


12 


6.28E- 


02 


6.29E- 


02 


5.65E- 


02 


1.41E- 


03 


6.94E - 


05 


1.67E- 


06 


5.46E- 


08 


23 


5.34E - 


03 


5.34E - 


03 


6.02E- 


03 


1.92E- 


03 


5.60E- 


05 


1.73E- 


06 


5.49E - 


08 


24 


2.20E- 


02 


2.20E- 


02 


2.26E- 


02 


2.79E- 


03 


6.28E- 


05 


1.68E- 


06 


5.47E - 


08 


47 


1.76E- 


03 


1.76E- 


03 


1.81E- 


03 


1.59E- 


03 


5.40E - 


05 


1.73E- 


06 


5.49E - 


08 


48 


7.78E - 


03 


7.78E - 


03 


7.83E - 


03 


8.18E- 


04 


5.07E- 


05 


1.68E- 


06 


5.47E - 


08 


95 


6.01E- 


04 


6.01E- 


04 


6.06E- 


04 


1.30E- 


03 


5.68E- 


05 


1.72E- 


06 


5.49E - 


08 


96 


2.75E - 


03 


2.75E- 


03 


2.75E- 


03 


3.19E- 


03 


6.43E - 


05 


1.66E- 


06 


5.47E - 


08 


N\k 







1 




10 




10 2 




10 3 




10 4 




10 5 




11 


1.10E- 


03 


1.12E- 


03 


2.36E- 


03 


1.59E- 


05 


3.03E- 


07 


3.58E- 


09 


4.59E- 


11 


12 


1.79E- 


03 


1.80E- 


03 


1.97E- 


03 


7.04E - 


05 


1.36E- 


06 


1.66E- 


08 


2.09E - 


10 


23 


1.54E- 


04 


1.55E- 


04 


2.07E- 


04 


3.87E- 


05 


1.95E- 


07 


1.97E- 


09 


2.56E- 


11 


24 


3.08E - 


04 


3.09E - 


04 


3.45E- 


04 


6.30E- 


05 


5.34E - 


07 


6.14E- 


09 


7.78E- 


11 


47 


2.44E - 


05 


2.44E - 


05 


2.64E- 


05 


2.15E- 


05 


7.43E - 


08 


7.59E- 


10 


1.05E- 


11 


48 


5.41E- 


05 


5.42E - 


05 


5.57E- 


05 


1.20E- 


05 


1.54E- 


07 


2.14E- 


09 


2.79E- 


11 


95 


4.07E - 


06 


4.08E - 


06 


4.16E- 


06 


1.75E- 


05 


1.16E- 


07 


2.00E - 


10 


3.91E- 


12 


96 


9.56E- 


06 


9.56E- 


06 


9.62E- 


06 


1.90E- 


05 


1.68E- 


07 


6.33E - 


10 


9.81E- 


12 



Table 10: Errors of the quadrature rule for integral ()5.2b[) with a = —1, (3 = 1/2 (top) and (3 = 
(bottom) 
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where C is independent of n. 

Unlike T n , U n is not uniformly bounded in n and x £ [—1,1]. However, 



\U n \\l,u = f \Un{x)\ 2 v 7 ! - x2 dx = jT sin 2 n0d0 = |. (A.5) 



On the other hand, 



f 1 T n (x) dx= T cos n9 sin Odd = i^ 1, lf n 18 even > ( A6 ) 
Jo I 0; otherwise. 

The trigonometric identity 

cos nO sin(m + 1)6* = ^ (sin(m + n + 1)0 + sin(m + 1 — n)6) 

implies 

f i (C/ m +„ + £7 m - n ) , if m > n - 1, 
T n *7 m = ^ (A. 7) 

In particular, we obtain for n > 1 



2xT^(x) = 2nT 1 (x)U n - 1 (x)=n[U n (x) + U n - 2 (x)], (A.8) 

il 



T n (z) = T n (x)[/o(x) = i[[/„(x)-[/ n _ 2 (x)]. (A.9) 



Finally, it holds 

T n {x)-T n (y) = 2 J2 Uj{x)Tn _ i _. {y) + Un _ i{x) = 2 J2u j (y)T n _ 1 _ j (x) + U n ^(y) (A.10) 
X V j=0 3=0 

which can be easily proven by induction on n. 
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